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Abstract 

Water structure, measured by the height of the first peak in oxygen-oxygen radial distributions, 
is convverged with respect to plane-wave basis energy cutoffs for ab initio molecular dynamics 
simulatinos, confirming the reliability of plane-wave methods. 



Given its ubiquitous presence both in our external and internal environments, and its 
active participation in numerous chemical processes, water is arguably the most important 
of all liquids. Yet despite extensive study and documentation water has been called the least 
understood of the known liquids with many of its properties considered anomalous.- 

In the thirty-five years since Franks' comprehensive assessment,- much progress has been 
made toward attaining a molecular-level understanding af water, particularly with regard to 
water structure.^ High accuracy scattering data collected over this time, including the most 
recent studies,- converge on a view that water is modestly structured. Accurate structural 
data enables quantitative assessment of computer simulations, which potentially can play a 
pivotal role in advancing our understanding not only of water and its anomalies, but also of 
solvation and interfacial water behavior.- We focus this study on ab initio molecular dynam- 
ics (AIMD) simulations of water and assess whether structural predictions are converged 
with respect to a particular choice of methodology. 

Most AIMD simulations of water-^'^'^'^'^'i^'ii'i^ using the BLYP^^ and PBE^^ functional 
report enhanced structure of liquid water relative to the most recent X-ray scattering data.- 
Over-structured water is predicted even though these simulations use disparate methods 
including varied pseudo-potentials, and simulation protocols. An exception occurs when 
the RPBE^ exchange-correlation functional is applied, for which the predicted water struc- 
ture nearly matches experiment.- All of these simulations cited share one common feature: 
expansion of the wavefunctions in a plane- wave basis set. 

A recent AIMD study by Lee and Tuckerman^^ (henceforth referred to as LT) apply an al- 
ternative basis set and report a less significant over-structuring of water relative to most prior 
simulations. Specifically, LT find that their simulations using the BLYP exchange-correlation 
functional^ yield an 0-0 pair correlation function g{r) below 3.0 at temperature 300 K, 
indicating a less significant over-structuring than earlier works .~'^'^'^'^'^>^>^>^ The authors 
employ a real-space method to expand the wavefunctions. Discrete Variable Representation 
(DVR), and stress the high degree of energy and force convergence achieved using that basis 
set. LT also speculated that effects of an incomplete plane wave basis set and consequent 
inadequate force convergence in previous woTk—^^^^^^^'^^^^ could in general be responsible 
for overstructured water g{r) published in the literature for multiple functionals, including 
not just BLYP but also PBE and PW91.i^ 



The discrepancy between DVR and existing plane wave predictions for BLYP water, 
attributed to under convergence of plane wave basis used in the literature/^ has motivated 
us to examine this issue for another functional widely used for water, namely PBE. We 
ask the question of what EC is required to converge water structure for AIMD simulations 
using plane wave basis set expansions. Toward this end, we present plane-wave based AIMD 
simulations of water (D2O), and systematically analyze convergence of the maximum height 
of the first peak in the oxygen-oxygen g{r) with respect to the wavefunction energy cutoff 
(EC). 

The simulations are made using the code VASP.— We study cells with 64 and 32 molecules 
at 1.0 g/cc density (assuming, for this purpose only, the proton mass) in the NVT (Nose- 
Hoover thermostat) ensemble at 300 K and 375 K, starting from several different initial 
configurations with F-point Brillouin zone sampling. We employ as a general model the 
PBE^^ exchange-correlation functional with plane-augmented waves applied to describe the 
interaction between core and valence electrons, a pseudo-potential which has been used 
extensively for many systems. The Born-Oppenheimer energy convergence criteria, time 
steps (0.25 or 0.5 ps) used, and total energy drifta^^ in these runs are listed in Table 1. 

Equilibration protocols used are these: run A started from the last configuration of a 
simulation using the empirical SPC/E^ force field (2 ps equilibration); B, from the end of 
A (0 ps); C, from the end of B (0 ps); D, from the end of A (6.2 ps); E, from a T=400 K 
simulation (12 ps); F, 8 ps into E (3 ps); G, end of E (0 ps); H, 8 ps into E (5 ps). 

Figure 1 illustrates equilibration of the thermostat for simulations conducted on 32 waters 
at temperatures of 375 K (trajectory B) and 300 K (trajectory D). (v) of each atom type 
(X =0 or H, for N =32 oxygens or 64 protons, using the deuterium mass for protons), 
scaled by Boltzmann's constant (/cb) and system temperature (T): 

X=N 

T= Y,im^x\vx\V{Sk^/2). (1) 
x=i 

Equilibration of both O and H atoms to the same temperature is achieved within 3- 
4 ps of simulation time using a single global Nose-Hoover thermostat. During the 12 ps 
equilibration of run E to ambient temperature from a high-temperature initial configuration 
(T=400 K), the g{r) peak value remained between 3.0 and 3.4 (not illustrated), in other 
words over-structured, until equilibrium was established. 



Figure 2a shows the 0-0 and H-H pair correlation functions for D2O at 375 K, with peak 
heights further hsted in Table 1. At this temperature, the pair-correlation functions do not 
become less overstructured when EC is increased by a factor of two beyond the VASP recom- 
mended energy cutoff of 400 eV for our oxygen and hydrogen PAW pseudopotentials. (Note 
that PAW simulations generally require a lower EC than norm-conserving pseudopotentials 
to yield converged water g{r).) 

At 300 K (see Fig. 2b) the slow dynamics of PBE water cause larger statistical errors 
in 0-0 peak height compared to the 375 K runs. The width and rise of the distribution 
functions are, in contrast, well converged. Faster dynamics in the 375 K runs result in more 
efficiently sampled simulations, making simulations at 375 K more relevant demonstrations 
of peak- height convergence with EC due to this statistical effect. Nevertheless, there is 
no reason to conclude that convergence behavior of AIMD simulations with respect to EC 
would be different at 375 K and at 300 K; the difference in thermal energy between the 
two temperatures is in fact negligible in comparison to the energy fluctuations present in an 
NVT simulation at these temperatures. 

Our goo{f) computed at T=300 K (Fig. 2b) are similar to those of Ref. [5] and 3, the 
former of which was performed at a slightly higher temperature of T=337 K. Considering the 
statistical uncertainties, it is more significant that at these temperatures the PBE functional 
consistently predicts considerable overstructuring of water. 

The case of EC=300 eV (Fig. 2b, Trajectory F in Table 1) is chosen to demonstrate 
the kind of deviations caused by a too low EC. For example, gooi'f') appears significantly 
less over-structured than converged, EC=400 or 500 eV simulations, and its first peak is 
shifted to a smaller r value at this EC. At least for the PBE functional, our observed trend 
does not support LT's speculation that raising EC to improve convergence may yield less 
overstructured water. 

A possible source of errors in AIMD simulations is the fidelity of the pseudo-potentials. 
We have verified that employing a harder set of potentials (PAW with nominal 700 eV cutoff) 
does not lead to changes in g{r). 

Finally, given the smooth and asymptotic character of convergence with respect to EC 



Fig. 1 in Ref. Il6l ). it is highly unlikely that going well beyond EC=1000 eV could result in 



any changes in structural properties. 



Our main conclusion is thus that it is crucial to assess convergence of the property of 
interest (here 5'oo('")) iii any kind of AIMD simulation. The property needs to be converged 
with respect to, for example, plane-wave energy cut-off, as well as fictitious mass, fc-point, 
type and hardness of pseudo-potentials, time-step, and system size. Our independent and 
well- converged calculations of the 0-0 pair correlation function for water, conducted while 
systematically varying EC and enforcing identical pseudopotentials, Brillouin zone sam- 
pling, and simulation protocol (Born-Oppenheimer as opposed to Car-Parrinello molecular 
dynamics), indicate that raising EC does not lessen overstructuring of PEE water at 1.0 
g/cc density. In contrast, convergence studies of different properties, like energy or pressure 
(see for example Ref. and I21I). may or may not result in a converged structure. 

To summarize, although the results of Lee and Tuckerman are interesting, and are sure to 
stimulate additional work on this important subject, we have shown that the overstructured 
water g{r)^s computed at T ~ 300 K using the PBE functional published in the literature 
are not artifacts due to underconverged plane-wave cutoffs. 
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FIG. 1: Atomic temperatures (T) in units Kelvin averaged for 2 ps time segments over 32 oxygen 
(O, depicted in black) and 64 hydrogen (X, red) atoms during the time course (t) in ps of a) 
Trajectory B, b) Trajectory D. The error bars represent one standard deviation in temperature for 
each atomic species at each time step. 
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FIG. 2: Pair-correlation functions g}m{r) (at shorter distances) and goo{f) (at longer distances) 
for D2O. (a) 32 molecule cell at 375 K: EC=400 (green), 800 (blue), and 1200 eV (black) (runs A, 
B, C in Table 1). (b) 64 molecules cell at 300 K: EC=300 (red), 500 (blue), 1000 (black) (runs F, 
G, H in Table 1). 
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TABLE I: Details of our AIMD/PBE simulations. Each trajectory {traj.) is identified by an al- 
phabetical label, the number of waters {Neater) in the simulation box, the wavefunction energy 
cutoff in eV {EC), the thermostat temperature in Kelvin {temp), the energy convergence tolerance 
between successive self-consistent iterations during electronic structure calculations {conv), simu- 
lation time in ps {t), timestep in ps {At), maximum value of the first peak in the oxygen-oxygen 
radial distribution function {goo{f)max.), and the drift in temperature per ps {drift). 



